“The 2020 Environmental Performance Index (EPI) provides a data-driven summary of the state of sustainability around the world. Using 32 performance indicators across 11 issue categories, the EPI ranks 180 countries on environmental health and ecosystem vitality. These indicators provide a gauge at a national scale of how close countries are to established environmental policy targets. The EPI offers a scorecard that highlights leaders and laggards in environmental performance and provides practical guidance for countries that aspire to move toward a sustainable future.” – 2020 Environmental Performance Index
Load used libraries.
library("reshape2")
library("ggplot2")
theme_set(theme_bw(base_size = 12) +
theme(rect = element_rect(fill = "transparent")))
library("leaflet")
library("rgdal")
library("htmlwidgets")
library("lavaan")
library("lavaanPlot")
Clear the workspace and set working directory.
rm(list = ls())
setwd("~/indblik/epi")
This helper function retrieves, reads and processes the EPI- and the GDP from the World Bank.
source("importEPI.R")
Run the helper function importEPI() and release the contents of the resulting list into the global environment.
epi.list <- importEPI()
## 'data.frame': 8234 obs. of 7 variables:
## $ code : int 4 24 8 784 32 51 28 36 40 31 ...
## $ iso : chr "AFG" "AGO" "ALB" "ARE" ...
## $ country : Factor w/ 179 levels "Afghanistan",..: 1 4 2 170 6 7 5 8 9 10 ...
## $ region : Factor w/ 8 levels "Asia-Pacific",..: 7 8 2 5 6 3 6 4 4 3 ...
## $ GDP2019 : num 1.93e+10 8.88e+10 1.53e+10 4.21e+11 4.45e+11 ...
## $ EPI.new : Factor w/ 46 levels "AGR","AIR","APE",..: 12 12 12 12 12 12 12 12 12 12 ...
## $ EPI.new.value: num 25.5 29.7 49 55.6 52.2 52.3 48.5 74.9 79.6 46.5 ...
list2env(epi.list, .GlobalEnv)
## <environment: R_GlobalEnv>
rm(epi.list)
Download the shape file with country borders and some affiliated data (area, population etc.) from thematicmapping, if it is absent. Note, this overwrites existing files!
if(!file.exists("2020/TM_WORLD_BORDERS-0.3.shp")){
print("File does not exist, downloading it now:")
download.file("thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip",
destfile = "world_shape_file.zip")
system("unzip -u world_shape_file.zip;
rm world_shape_file.zip Readme.txt")
}
Then, read this shape file into a SpatialPolygonsDataFrame and inspect it.
world.spdf <- readOGR(dsn = "2020", layer = "TM_WORLD_BORDERS-0.3",
verbose = TRUE)
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/r/indblik/epi/2020", layer: "TM_WORLD_BORDERS-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings: POP2005
summary(world.spdf)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -180 180.0000
## y -90 83.6236
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
## FIPS ISO2 ISO3 UN
## Length:246 Length:246 Length:246 Min. : 4.0
## Class :character Class :character Class :character 1st Qu.:215.0
## Mode :character Mode :character Mode :character Median :429.0
## Mean :431.8
## 3rd Qu.:650.5
## Max. :894.0
## NAME AREA POP2005 REGION
## Length:246 Min. : 0.0 Length:246 Min. : 0.00
## Class :character 1st Qu.: 44.5 Class :character 1st Qu.: 2.00
## Mode :character Median : 5515.5 Mode :character Median : 19.00
## Mean : 52696.1 Mean : 65.43
## 3rd Qu.: 34708.8 3rd Qu.:142.00
## Max. :1638094.0 Max. :150.00
## SUBREGION LON LAT
## Min. : 0.00 Min. :-178.13 Min. :-80.4460
## 1st Qu.: 14.00 1st Qu.: -50.16 1st Qu.: -0.3025
## Median : 30.00 Median : 17.66 Median : 16.5110
## Mean : 54.84 Mean : 13.29 Mean : 16.4075
## 3rd Qu.: 61.00 3rd Qu.: 50.01 3rd Qu.: 39.1067
## Max. :155.00 Max. : 179.22 Max. : 78.8300
Unfortunately, as soon as it comes to map the world, geopolitics come to the fore. This is reflected in the mismatches between country lists due to the different recognition statuses of these countries. Luckily, an ISO standard exists that allocates a unique 3-letter code to every country.
For one, there are distinctly more countries in the shape file:
length(world.spdf$NAME)
## [1] 246
length(epi$country)
## [1] 179
These mismatches stem mainly from different naming schemes…
sort(setdiff(world.spdf$NAME, epi$country))
## [1] "Åland Islands"
## [2] "American Samoa"
## [3] "Andorra"
## [4] "Anguilla"
## [5] "Antarctica"
## [6] "Aruba"
## [7] "Bermuda"
## [8] "Bouvet Island"
## [9] "British Indian Ocean Territory"
## [10] "British Virgin Islands"
## [11] "Burma"
## [12] "Cape Verde"
## [13] "Cayman Islands"
## [14] "Christmas Island"
## [15] "Cocos (Keeling) Islands"
## [16] "Congo"
## [17] "Cook Islands"
## [18] "Democratic Republic of the Congo"
## [19] "Falkland Islands (Malvinas)"
## [20] "Faroe Islands"
## [21] "French Guiana"
## [22] "French Polynesia"
## [23] "French Southern and Antarctic Lands"
## [24] "Gibraltar"
## [25] "Greenland"
## [26] "Guadeloupe"
## [27] "Guam"
## [28] "Guernsey"
## [29] "Heard Island and McDonald Islands"
## [30] "Holy See (Vatican City)"
## [31] "Hong Kong"
## [32] "Iran (Islamic Republic of)"
## [33] "Isle of Man"
## [34] "Jersey"
## [35] "Korea, Democratic People's Republic of"
## [36] "Korea, Republic of"
## [37] "Lao People's Democratic Republic"
## [38] "Libyan Arab Jamahiriya"
## [39] "Liechtenstein"
## [40] "Macau"
## [41] "Martinique"
## [42] "Mayotte"
## [43] "Micronesia, Federated States of"
## [44] "Monaco"
## [45] "Montserrat"
## [46] "Nauru"
## [47] "Netherlands Antilles"
## [48] "New Caledonia"
## [49] "Niue"
## [50] "Norfolk Island"
## [51] "Northern Mariana Islands"
## [52] "Palau"
## [53] "Palestine"
## [54] "Pitcairn Islands"
## [55] "Puerto Rico"
## [56] "Republic of Moldova"
## [57] "Reunion"
## [58] "Saint Barthelemy"
## [59] "Saint Helena"
## [60] "Saint Kitts and Nevis"
## [61] "Saint Martin"
## [62] "Saint Pierre and Miquelon"
## [63] "San Marino"
## [64] "Somalia"
## [65] "South Georgia South Sandwich Islands"
## [66] "Svalbard"
## [67] "Swaziland"
## [68] "Syrian Arab Republic"
## [69] "Taiwan"
## [70] "The former Yugoslav Republic of Macedonia"
## [71] "Tokelau"
## [72] "Turks and Caicos Islands"
## [73] "Tuvalu"
## [74] "United Republic of Tanzania"
## [75] "United States"
## [76] "United States Minor Outlying Islands"
## [77] "United States Virgin Islands"
## [78] "Wallis and Futuna Islands"
## [79] "Western Sahara"
## [80] "Yemen"
… but also from the fact that the EPI is missing for some countries. Countries without an EPI are mostly miniature states (e.g. Monaco, Liechtenstein or Andorra), states currently in a(n) (armed) conflict (Syria, Palestine or Taiwan) or several islands states. Some of these islands actually belong to a country; e.g., Greenland and Svalbard that belong to Denkmark Norway, respectively. For now, keep it as is although it would be reasonable to merge these islands with their mainlands.
sort(setdiff(epi$country, world.spdf$NAME))
## [1] "Cabo Verde" "Dem. Rep. Congo"
## [3] "Eswatini" "Iran"
## [5] "Laos" "Micronesia"
## [7] "Moldova" "Myanmar"
## [9] "North Macedonia" "Republic of Congo"
## [11] "South Korea" "Tanzania"
## [13] "United States of America"
Now, merge the epi data.frame and the world.spdf SpatialPolygonsDataFrame by the country ISO-3 so that the latter holds all relevant information to render an interactive choropleth map.
world.spdf <- merge(world.spdf,
epi[, (colnames(epi) %in% c("iso", "EPI.new", "region",
"GDP2019"))],
by.x = "ISO3", by.y = "iso")
names(world.spdf)
## [1] "ISO3" "FIPS" "ISO2" "UN" "NAME" "AREA"
## [7] "POP2005" "REGION" "SUBREGION" "LON" "LAT" "EPI.new"
## [13] "region" "GDP2019"
For readability, population numbers are encoded as millions (M) and for meaningfulness, as numbers instead of characters.
world.spdf@data$POP2005[which(world.spdf@data$POP2005 == 0)] <- NA
world.spdf@data$POP2005 <- as.integer(world.spdf@data$POP2005) / 1000000 %>%
round(2)
Now, let’s render an interactive map with R leaflet to display the EPI around the globe.
First, create a gradient color scheme for the EPI…
epi.cols <- colorNumeric(palette = "viridis", domain = world.spdf@data$EPI.new,
na.color = "gray", reverse = TRUE)
…and then a discrete color scheme for the regions.
region.cols <- colorFactor(palette = "plasma", domain = world.spdf@data$region,
na.color = "gray")
Assemble data for tooltips that emerge when hovering over the html-ified map.
tooltips <- paste("Country:", world.spdf@data$NAME, "<br/>",
"Region:", world.spdf@data$region, "<br/r>",
"GDP:", round(world.spdf@data$GDP2019 / 1000000000,
digits = 1), "B US$", "<br/r>",
"Population:", round(world.spdf@data$POP2005, 2), "M",
"<br/>",
"EPI:", world.spdf@data$EPI.new) %>%
lapply(htmltools::HTML)
Finally, use the previous steps to generate a leaflet choropleth map.
epi.map <- leaflet(world.spdf) %>%
addTiles() %>%
setView(lat = 20, lng = 0, zoom = 1.5) %>%
addPolygons(fillColor = ~ epi.cols(EPI.new), stroke = TRUE, group = "EPI",
fillOpacity = 0.9, color = "white", weight = 0.5,
label = tooltips,
labelOptions = labelOptions(style = list("font-weight" = "normal",
padding = "3px 8px"),
textsize = "13px",
direction = "auto")) %>%
addPolygons(fillColor = ~ region.cols(region), stroke = TRUE,
group = "Region", fillOpacity = 0.5, color = "white",
weight = 0.5, label = tooltips,
labelOptions = labelOptions(style = list("font-weight" = "normal",
padding = "3px 8px"),
textsize = "13px",
direction = "auto")) %>%
addLegend(pal = epi.cols, values = ~ EPI.new, opacity = 0.9,
title = "EPI", position = "topleft", group = "EPI") %>%
addLegend(pal = region.cols, values = ~ region, opacity = 0.5,
title = "", position = "topleft", group = "Region") %>%
addLayersControl(baseGroups = c("EPI", "Region"),
options = layersControlOptions(collapsed = FALSE)) %>%
hideGroup("Region")
epi.map